//Redfish finite-state, continuous-time, age-structured model
//Tim Miller January 2008
DATA_SECTION
//READ DATA FROM INPUT FILE "FSCT_GBYT.DAT"
//NO SPACES in (,)
  init_number myPI
  init_int styr
  init_int endyr
  init_int n_N0s
  init_int n_Fs
  init_int n_cohorts
  //need to add these indicators of where to start in catchdate/catchlength data
  init_ivector N0_indicators(1,n_cohorts)
  //  init_ivector selectivity_types(1,3) //1: survey 1, 2: survey 2, 3: commercial. If type = 0: selectivity at age, type=1: logistic
  
  init_ivector cohorts(1,n_cohorts)
  init_vector cohort_jbdates(1,n_cohorts)
  init_ivector cohort_start_periods(1,n_cohorts) //index period to start for each cohort
  init_vector cohort_ages_at_start_periods(1,n_cohorts) //ages at the beginning of the start periods
  init_ivector survey_cohort_data_starts(1,n_cohorts) //index in aged survey data to start for each cohort
  init_ivector catch_cohort_data_starts(1,n_cohorts) //index in aged catch data to start for each cohort
  
  init_int nperiods
  init_int ncruises
  init_int nlengths
  init_int stlen
  init_int endlen
  init_ivector len_cl(1,nlengths)
  vector lengths(1,nlengths)
  !! lengths = len_cl + 0.5; // Assume len_cl index lower bound of 1 cm interval...
  init_int n_survey_agedata
  init_ivector n_survey_by_cohort(1,n_cohorts)
  init_int n_catch_agedata
  init_ivector n_catch_by_cohort(1,n_cohorts)
  
  //ageing probabilities for each survey...
  init_matrix survey_natlength(1,ncruises,stlen,endlen)
  init_matrix survey_nagedatlength(1,ncruises,stlen,endlen)
  
  //ageing probabilities for landings are in terms of years...
  init_matrix catch_natlength(styr,endyr,stlen,endlen) 
  init_matrix catch_patlength(styr,endyr,stlen,endlen) 
  init_matrix catch_nagedatlength(styr,endyr,stlen,endlen)
  
  //periodinfo_idata has 5 rows of period-specific (columns) integer data: year, month, survey season, cruise indicator, Findicator
  init_int nperiodinfo_idata
  init_imatrix periodinfo_idata(1,nperiodinfo_idata,1,nperiods)
  //periodinfo_ddata has 8 rows of period-specific (columns) double data: landmtyr, catch_nlengthed, survey_nlengthed, survey_ntotal,
  //towsperyr, startj, timeyr, yearlengths
  init_int nperiodinfo_ddata
  init_matrix periodinfo_ddata(1,nperiodinfo_ddata,1,nperiods)
  
  //these are numbers aged by cohort (rows) and period (columns)
  init_imatrix periodinfo_n_survey_aged(1,n_cohorts,1,nperiods)
  init_imatrix periodinfo_n_catch_aged(1,n_cohorts,1,nperiods)
  
  // THESE NEED TO BE in chronological order of catchdate
  init_vector survey_jcatchdate(1,n_survey_agedata)
  init_vector survey_catchlength(1,n_survey_agedata)

//THESE NEED TO BE in chronological order of catchdate
  init_vector catch_catchlength(1,n_catch_agedata)

  vector period_ddata(1,8);
  ivector period_idata(1,5);
	
INITIALIZATION_SECTION
 //PROVIDE INITIAL PARAMETER VALUES
 //natural mortality
 log_M -1.609438
 //log_M -2.995732
 
 //log(recruitment)
 log_N0 7.0

 log_qs -9.0
 
 //survey parameters
 survey1_selectivity_pars -0.69
 survey2_selectivity_pars -0.69
 catch_selectivity_pars -0.69

 //commercial catch parameters
 log_F -1.6

 //growth parameters
 log_Linf 3.5971872
 log_k_LVB -1.8004082
 a0 0.2435075
 log_sig_LVB -2.0652059

 //lenwt parameters
 log_alpha -11.78885
 log_beta 1.17088
 log_sig_lw -2.05294

PARAMETER_SECTION
 //DECLARE MODEL PARAMETERS AND VARIABLES
 init_bounded_number log_M(-10.,10.,-1)
 // init_bounded_vector log_N0(1,n_N0s,-35.0,50.5,1)
 init_vector log_N0(1,n_N0s,1)

 // init_bounded_vector log_qs(1,2,-35.0,50.0,2)
 // init_bounded_vector log_qs(1,1,-25.0,-5.0,2)
 init_vector log_qs(1,1,4)
 
 // init_vector survey1_selectivity_pars(1,2,4)
 // init_vector survey2_selectivity_pars(1,1,4)
 init_vector survey1_selectivity_pars(1,1,-1)
 init_vector survey2_selectivity_pars(1,1,-1)
 init_vector catch_selectivity_pars(1,1,-1)
 // init_vector catch_selectivity_pars(1,1,4)

 // init_bounded_vector log_F(1,n_Fs,-35.0,50.5,3)
 init_vector log_F(1,n_Fs,3)

 init_number log_Linf(-1)
 init_number log_k_LVB(-1)
 init_number a0(-1)
 init_number log_sig_LVB(-1)

 init_number log_alpha(-1)
 init_number log_beta(-1)
 init_number log_sig_lw(-1)

 number M
 vector N0(1,n_N0s)
 number current_N0
 vector qs(1,2)
 vector F(1,n_Fs)
 number Linf
 number k_LVB
 number sig_LVB
 number alpha
 number beta
 number sig_lw

 number ll_catch_age
 number p_for_catch_data
 number ll_survey_age
 number ll_notaged
 number ll_first
 vector temp_ddat(1,7)
 matrix nhats(1,n_cohorts,1,2)
 vector wt_at_len_hat(1,nlengths);
		 
 objective_function_value total_ll

RUNTIME_SECTION
 //no tabs!!
 convergence_criteria 1e-6;

PRELIMINARY_CALCS_SECTION
 //CHECK INPUT DATA
 //must be space(s) at the beginning of each line??!!!!
 cout << "PI: " << PI << endl;
 cout << "myPI: " << myPI << endl;
 cout << "START YEAR: "<<styr<< endl;
 cout << "END YEAR: "<<endyr<< endl;
 cout << "n cohorts: "<<n_cohorts<<endl;
 cout << "n N0 pararameters: "<<n_N0s<< endl;
 cout << "n F pararameters: "<<n_Fs<< endl;
 cout << "cohort years: "<< cohorts << endl;
 cout << "n periods (period info): "<< nperiods<< endl;
 cout << "start length: "<< stlen << endl;
 cout << "end length: "<< endlen << endl;
 cout << "nlength classes: "<< nlengths << endl;
// cout << "N0 indicators: "<< N0_indicators << endl;
 cout << "n survey age data: "<< n_survey_agedata<< endl;
 cout << "n_catch_agedata: " << n_catch_agedata << endl;
// cout << "catch age data length at catch: "<< catch_catchlength << endl;
 beta = mfexp(log_beta);
 sig_lw = mfexp(log_sig_lw);
 for(int j = 1; j <= nlengths; j++)
 {
   wt_at_len_hat(j) = alpha * pow(lengths(j),beta);
 }
 
TOP_OF_MAIN_SECTION
 //ALLOCATE SPACE IN READ-WRITE MEMORY
 arrmblsize=200004000;
 gradient_structure::set_GRADSTACK_BUFFER_SIZE(2004000);
 gradient_structure::set_CMPDIF_BUFFER_SIZE(60000000);
 gradient_structure::set_NUM_DEPENDENT_VARIABLES(20000);

GLOBALS_SECTION

PROCEDURE_SECTION
  Initialize_Model();
  for(int i = 1;i <= n_cohorts; i++)
  {
	  total_ll += get_ll_cohort(i);
	  if(last_phase())
    {
		   nhats(i) = get_nhat_cohort(i);
	  }
  }

FUNCTION Initialize_Model
 //NOTE: N0 is recruitment (numbers at birth) less those ending up aged in either survey or landings
 M = mfexp(log_M);
 N0 = mfexp(log_N0);
 F = mfexp(log_F);
 qs(1) = mfexp(log_qs(1));
 qs(2) = 0.0;
 //qs(2) = mfexp(log_qs(2));
 Linf = mfexp(log_Linf);
 k_LVB = mfexp(log_k_LVB);
 sig_LVB = mfexp(log_sig_LVB);
 alpha = mfexp(log_alpha);
 beta = mfexp(log_beta);
 sig_lw = mfexp(log_sig_lw);

 total_ll = 0.0;
 ll_catch_age = 0.0;
 ll_survey_age = 0.0;

 if(active(log_beta))
 {
   for(int j = 1; j <= nlengths; j++)
   {
 	   wt_at_len_hat(j) = alpha * pow(lengths(j),beta);
   }
 }

FUNCTION dvariable get_ll_cohort(int& cohort_ind) 
  RETURN_ARRAYS_INCREMENT();
	/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	//Returns variable "ll_cohort"
		
	//periodinfo_ddata matrix rows: period_landmtyr, periodinfo_catch_nlengthed, periodinfo_survey_nlengthed, periodinfo_survey_ntotal,
	//periodinfo_towsperyr, periodinfo_startj, periodinfo_timeyr, periodinfo_yearlengths
	//
	//periodinfo_idata matrix rows: periodinfo_year, periodinfo_month, periodinfo_survey_season, periodinfo_cruise_indicator, 
	//periodinfo_Findicators
	
	//NOTE!!!: catch_length, survey_jcatchdate, survey_length must be in chronological (catch date) order
	//NOTE!!!: cohort_start_ind needs to be the interval where the cohort birthdate occurs!!!!
	
	dvariable ll_catch = 0.0, ll_survey =0.0, ll_notaged = 0.0, ll_cohort = 0.0;
	dvariable mean_wt_hat = 0.0, catch_Nhat = 0.0, F_catch = 0.0;
	dvar_vector period_probs(1,4);
	dvar_vector surv_pars(1,8);
	dvar_vector cat_pars(1,4);
	int n_catch, n_survey, current_length_class = 0;
	dvariable cum_p_surv = 1.0, current_length = 0.0, cum_page_s1 = 0.0, cum_page_s2 = 0.0, cum_page_catch = 0.0;
	double time = 0.0, current_age = 0.0;
	double current_fish_length = 0.0, current_fish_jcdate = 0.0;
	//starting index for vector of fish indices (length is # of all fish in all cohorts)
	int current_fish_ind_catch = catch_cohort_data_starts(cohort_ind);
	int current_fish_ind_survey = survey_cohort_data_starts(cohort_ind);
	
	for(int i=cohort_start_periods(cohort_ind);i<=nperiods;i++){
		//period specific info 5 idata and 8 ddata
		period_idata = column(periodinfo_idata, i);
		period_ddata = column(periodinfo_ddata, i);
		if(i==cohort_start_periods(cohort_ind)){
			if(cohort_ages_at_start_periods(cohort_ind) < 0){
				current_age = 0.0;
				//need to shorten the first interval for the birthdate of the cohort
				period_ddata(7) = period_ddata(7) + (period_ddata(6) - cohort_jbdates(cohort_ind))/period_ddata(8);
				//then change the start date
				period_ddata(6) = cohort_jbdates(cohort_ind);
			}
			else current_age = cohort_ages_at_start_periods(cohort_ind); //cohorts starting before first year of model
		}
		F_catch = F(period_idata(5));

		//number of fish caught in current period for this cohort
		n_catch = periodinfo_n_catch_aged(cohort_ind,i);
		n_survey = periodinfo_n_survey_aged(cohort_ind,i);		
		
		mean_wt_hat = sum(elem_prod(wt_at_len_hat,extract_row(catch_patlength,periodinfo_idata(1,i))));
		catch_Nhat = 1000.0 * periodinfo_ddata(1,i)/mean_wt_hat;
					
		if(n_catch> 0){			
			//get lengths to use for current period
			for(int j = 1; j <= n_catch; j++){
				current_fish_length = catch_catchlength(current_fish_ind_catch);
				ll_catch += get_ll_catch_fish(period_ddata(7),F_catch,cum_p_surv,current_fish_length,
					catch_Nhat,catch_natlength,catch_nagedatlength,
					survey_natlength,survey_nagedatlength, current_age);
				
				current_fish_ind_catch += 1;
			}
		}
	// if(i==89) cout << "i: " << i << " here 3" << endl;
		if(n_survey> 0){
			//get lengths and catchdates to use for current period
			for(int j = 1; j <= n_survey; j++){
				current_fish_length = survey_catchlength(current_fish_ind_survey);
				current_fish_jcdate = survey_jcatchdate(current_fish_ind_survey);
			
				ll_survey += get_ll_survey_fish(F_catch, cum_p_surv, current_fish_length, current_fish_jcdate, 
					catch_Nhat,catch_natlength,catch_nagedatlength,
					survey_natlength,survey_nagedatlength, current_age);
				
				if(current_fish_ind_survey ==1){
					ll_first = ll_survey;
					temp_ddat(1) = current_fish_length;
					temp_ddat(2) = current_fish_jcdate;
					temp_ddat(3) = qs(1);
					temp_ddat(4) = qs(2);
					temp_ddat(5) = F_catch;
					temp_ddat(6) = cum_p_surv;
					temp_ddat(7) = catch_Nhat;
				}
				current_fish_ind_survey += 1;
			}
		}
	// if(i==89) cout << "i: " << i << " here 4" << endl;
		
		current_length = Linf * (1.0 - mfexp(-k_LVB * (current_age - a0)));		
		if(current_length >= 0.5){
			for(int j = 1; j <= nlengths; j++){
				if((current_length>= double(lengths(j)) -0.5) & (current_length < double(lengths(j)) + 0.5)){
					current_length_class = lengths(j);
					break;
				}
			}
		}
		
		//not first period
		time = period_ddata(7);
		cat_pars  = get_cat_pars(period_idata(1),current_length_class,catch_Nhat,period_ddata(2),F_catch,
			                      catch_natlength,catch_nagedatlength);
		surv_pars = get_surv_pars(period_idata(3),current_length_class,period_idata(4), period_ddata(5), 
                    period_ddata(4), period_ddata(3), qs, survey_natlength, survey_nagedatlength);

		int survey1_selectivity_type = 0, survey2_selectivity_type = 0, catch_selectivity_type = 0;
		surv_pars(1) = surv_pars(1) * get_s1_selectivity_age(survey1_selectivity_type, survey1_selectivity_pars, current_age);
		surv_pars(5) = surv_pars(5) * get_s2_selectivity_age(survey2_selectivity_type, survey2_selectivity_pars, current_age);
		cat_pars(1)  = cat_pars(1) * get_catch_selectivity_age(catch_selectivity_type, catch_selectivity_pars, current_age);
		period_probs = get_probs(time,surv_pars,cat_pars);
			
		cum_page_s1 += cum_p_surv * period_probs(2);
		cum_page_s2 += cum_p_surv * period_probs(3);
		
		cum_page_catch += cum_p_surv * period_probs(4);
		cum_p_surv *= period_probs(1);
		current_age += time;
	}
	ll_notaged = N0(N0_indicators(cohort_ind)) * log(1.0 - cum_page_s1 - cum_page_s2 - cum_page_catch);
	ll_cohort = ll_catch + ll_survey + ll_notaged;
	ll_cohort += gammln(N0(N0_indicators(cohort_ind)) + double(n_catch_by_cohort(cohort_ind)) + double(n_survey_by_cohort(cohort_ind)) + 1.0);
	ll_cohort -= gammln(N0(N0_indicators(cohort_ind)) + 1.0);
	ll_cohort -= gammln(double(n_catch_by_cohort(cohort_ind))+ 1.0);
	ll_cohort -= gammln(double(n_survey_by_cohort(cohort_ind)) + 1.0);
	ll_cohort *= -1.0;
	
		// cout << "ll_catch: " << ll_catch << endl;
		// cout << "ll_survey: " << ll_survey << endl;
		// cout << "ll_notaged: " << ll_notaged << endl;
		// cout << "cum_page_s1: " << cum_page_s1 << endl;
		// cout << "cum_page_s2: " << cum_page_s2 << endl;
		// cout << "cum_page_catch: " << cum_page_catch << endl;
		// cout << "surv_pars: " << surv_pars << endl;
		// cout << "qs: " << qs << endl;
		
		// cout << "N0_unobs: " << N0(N0_indicators(cohort_ind)) << endl;
//		cout << "N0_catch: " << double(n_catch_by_cohort(cohort_ind)) << endl;
//		cout << "N0_survey: " << double(n_survey_by_cohort(cohort_ind)) << endl;
//		cout << "log_fact_N0_all: " << gammln(N0(N0_indicators(cohort_ind)) + double(n_catch_by_cohort(cohort_ind)) + double(n_survey_by_cohort(cohort_ind)) + 1.0) << endl;
		// cout << "log_fact_N_unobs: " << gammln(N0(N0_indicators(cohort_ind)) + 1.0) << endl;
		// cout << "log_fact_N_catch: " << gammln(double(n_catch_by_cohort(cohort_ind))+ 1.0) << endl;
		// cout << "log_fact_N_survey: " << gammln(double(n_survey_by_cohort(cohort_ind)) + 1.0) << endl;
		// cout << "diff in log_N: " << gammln(N0(N0_indicators(cohort_ind)) + double(n_catch_by_cohort(cohort_ind)) + double(n_survey_by_cohort(cohort_ind)) + 1.0) - gammln(N0(N0_indicators(cohort_ind)) + 1.0) - gammln(double(n_catch_by_cohort(cohort_ind))+ 1.0) - gammln(double(n_survey_by_cohort(cohort_ind)) + 1.0) << endl;
	RETURN_ARRAYS_DECREMENT();
	return ll_cohort;
	
// FUNCTION dvariable get_ll_catch_fish(double& time, dvar_vector& qs, dvariable& m, dvariable& 
//                                                            F_catch, dvariable& cum_p_surv, double& catchlength, int& nlengths, ivector& lengths, dvector& period_ddata, ivector& period_idata, dvariable& N_hat, dmatrix& cat_natlength,dmatrix& cat_nagedatlength, dmatrix& surv_natlength,dmatrix& surv_nagedatlength, double& current_age)
				//         get_ll_catch_fish(period_ddata(7),F_catch,                     cum_p_surv,current_fish_length,nlengths,lengths,
FUNCTION dvariable get_ll_catch_fish(double& time, dvariable& F_catch, dvariable& cum_p_surv, 
                  double& catchlength, dvariable& N_hat, dmatrix& cat_natlength,dmatrix& cat_nagedatlength, 
                  dmatrix& surv_natlength,dmatrix& surv_nagedatlength, double& current_age)
  RETURN_ARRAYS_INCREMENT();
	/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	//Returns variable "ll_catch_fish"
	//
	//period_ddata: landmtyr, catch_nlengthed, survey_nlengthed, survey_ntotal, towsperyr, startj, timeyr, yearlengths
	//period_idata matrix rows: year, month, survey_season, cruise_indicator, Findicators
	//
	//Requires: duration of period (time), survey q's (qs), natural mortality (m), commericial F (F_catch),
	//cumulative survival probabiltiy up to beginning of period (cum_p_surv), 
	//
	//lengths of fish at capture (cat_catchlengths), number of length classes (nlengths),
	//length classes (lengths), current catch Nhat (N_hat)
	//numbers at length matrix for catch (cat_natlength),numbers aged at length for catch (cat_nagedatlength),
	//numbers at length matrix for surveys (surv_natlength),
	//numbers aged at length matrix for surveys (surv_nagedatlength)
	/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

	dvariable ll_catch_fish = 0.0;
	dvar_vector surv_pars(1,8);
	dvar_vector cat_pars(1,4);
	dvar_vector period_probs(1,4);
	double len, cat_nlengthed, surv_ntot, surv_nlengthed, towsperyr;
	int tmp_len_cl, yr, surv_seas, crs_ind;
	
	cat_nlengthed = period_ddata(2);
	surv_nlengthed = period_ddata(3);
	surv_ntot = period_ddata(4);
	towsperyr = period_ddata(5);
	
	yr = period_idata(1);
	surv_seas = period_idata(3);
	crs_ind = period_idata(4);;
	len = catchlength;
	for(int j = 1; j <= nlengths; j++){
		if((len >= double(lengths(j)) - 0.5) & (len < double(lengths(j)) + 0.5)){
			len_cl = lengths(j);
			break;
		}
	}
//	cout << "Here catch 1" << endl;
	surv_pars = get_surv_pars(surv_seas,len_cl,crs_ind,towsperyr,surv_ntot,surv_nlengthed,qs,surv_natlength,surv_nagedatlength);
	cat_pars = get_cat_pars(yr,len_cl,N_hat,cat_nlengthed,F_catch,cat_natlength,cat_nagedatlength);
	int survey1_selectivity_type = 0, survey2_selectivity_type = 0, catch_selectivity_type = 0;
//	cout << "Here catch 2" << endl;
	surv_pars(1) = surv_pars(1) * get_s1_selectivity_age(survey1_selectivity_type, survey1_selectivity_pars, current_age);
	surv_pars(5) = surv_pars(5) * get_s2_selectivity_age(survey2_selectivity_type, survey2_selectivity_pars, current_age);
	cat_pars(1) = cat_pars(1) * get_catch_selectivity_age(catch_selectivity_type, catch_selectivity_pars, current_age);
	period_probs = get_probs(time,surv_pars,cat_pars);
	
	ll_catch_fish = log(cum_p_surv) + log(period_probs(4));
//	cout << "Here catch 4" << endl;
	RETURN_ARRAYS_DECREMENT();
	return ll_catch_fish;

	
FUNCTION dvariable get_ll_survey_fish(dvariable& F_catch, dvariable& cum_p_surv,double& catchlength,
                   double& jcatchdate, dvariable& N_hat, dmatrix& cat_natlength,dmatrix& cat_nagedatlength,dmatrix& surv_natlength,dmatrix& surv_nagedatlength, double& current_age)
  RETURN_ARRAYS_INCREMENT();
	/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	//Returns variable "ll_survey_fish"
	//
	//period_ddata: landmtyr, catch_nlengthed, survey_nlengthed, survey_ntotal, towsperyr, startj, timeyr, yearlengths
	//period_idata matrix rows: year, month, survey_season, cruise_indicator, Findicators
	//
	//Requires: natural mortality (m), 
	//cumulative survival probabiltiy up to beginning of period (cum_p_surv), 
	//length at capture (catchlength), catch date (jcatchdate),
	//number of length classes (nlengths), length classes (lengths), start of period (startj), length of current year (yrlen), current year (yr), 
	//current catch Nhat (N_hat), number lengthed in catch (cat_nlengthed) , commericial F (F_catch),
	//numbers at length matrix for catch (cat_natlength),numbers aged at length for catch (cat_nagedatlength),
	//current survey (surv_seas), current cruise (crs_ind), towsperyr, total number caught in survey (surv_ntot), 
	//number lengthed in survey (surv_nlengthed), survey q's (qs), numers at length matrix for surveys (surv_natlength),
	//numbers aged at length matrix for surveys (surv_nagedatlength)
	/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

	dvariable ll_survey_fish = 0.0;
	dvar_vector surv_pars(1,8);
	dvar_vector cat_pars(1,4);
	dvar_vector period_probs(1,4);
	int len_cl, yr, surv_seas, crs_ind;
	double len, cat_nlengthed, surv_ntot, surv_nlengthed, towsperyr, startj, yrlen, time;
	
	cat_nlengthed = period_ddata(2);
	surv_nlengthed = period_ddata(3);
	surv_ntot = period_ddata(4);
	towsperyr = period_ddata(5);
	startj = period_ddata(6);
	yrlen = period_ddata(8);
	
	yr = period_idata(1);
	surv_seas = period_idata(3);
	crs_ind = period_idata(4);;
	
	len = catchlength;
	for(int j = 1; j <= nlengths; j++){
		if((len >= double(lengths(j)) - 0.5) & (len < double(lengths(j)) + 0.5)){
			len_cl = lengths(j);
			break;
		}
	}
			
	time = (jcatchdate - startj)/yrlen;
			
	surv_pars = get_surv_pars(surv_seas,len_cl,crs_ind,towsperyr,surv_ntot,surv_nlengthed,qs,surv_natlength,surv_nagedatlength);
	cat_pars = get_cat_pars(yr,len_cl,N_hat,cat_nlengthed,F_catch,cat_natlength,cat_nagedatlength);
	int survey1_selectivity_type = 0, survey2_selectivity_type = 0, catch_selectivity_type = 0;
	surv_pars(1) = surv_pars(1) * get_s1_selectivity_age(survey1_selectivity_type, survey1_selectivity_pars, current_age);
	surv_pars(5) = surv_pars(5) * get_s2_selectivity_age(survey2_selectivity_type, survey2_selectivity_pars, current_age);
	cat_pars(1) = cat_pars(1) * get_catch_selectivity_age(catch_selectivity_type, catch_selectivity_pars, current_age);
	period_probs = get_probs(time,surv_pars,cat_pars);
			
	ll_survey_fish = log(cum_p_surv) + log(period_probs(1));
	if(surv_seas == 1){
		ll_survey_fish += log(surv_pars(1)) + log(surv_pars(2)) + log(surv_pars(3)) + log(surv_pars(4));
	}
	if(surv_seas == 2){
		ll_survey_fish += log(surv_pars(5)) + log(surv_pars(6)) + log(surv_pars(7)) + log(surv_pars(8));
	}
	RETURN_ARRAYS_DECREMENT();
	return ll_survey_fish;

FUNCTION dvar_vector get_probs(double& time, dvar_vector& surv_psamp, 
               dvar_vector& cat_psamp)
  RETURN_ARRAYS_INCREMENT();
	/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	//Returns vector "probs": survived, aged survey 1, aged survey 2, aged catch at end of "time"
	//
	//Requires: time, survey parameters (surv_psamp), catch parameters (cat_psamp), natural mortality (m)
	/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	
	dvar_vector probs(1,4);
	dvariable aged1, aged2, aged_catch, hazard, multiplier;
	
	aged1 = surv_psamp(1) * surv_psamp(2) * surv_psamp(3) * surv_psamp(4);
	aged2 = surv_psamp(5) * surv_psamp(6) * surv_psamp(7) * surv_psamp(8);
	aged_catch = cat_psamp(1) * cat_psamp(2) * cat_psamp(3) * cat_psamp(4);
	
	hazard = surv_psamp(1) + surv_psamp(5) + cat_psamp(1) + M;
	probs(1) = mfexp(-hazard * time);
	multiplier = (1.0 - probs(1))/ hazard;
	probs(2) = aged1 * multiplier;
	probs(3) = aged2 * multiplier;
	probs(4) = aged_catch * multiplier;
	
	RETURN_ARRAYS_DECREMENT();
	return probs;

FUNCTION dvar_vector get_surv_pars(int& surv_seas,int& len_cl,int& crs_ind,double& towsperyr,double& surv_ntot,double& surv_nlengthed, dvar_vector& qs, dmatrix& natlength, dmatrix& nagedatlength)
// FUNCTION dvar_vector get_surv_pars(int& surv_seas,int& len_cl,int& crs_ind,double& towsperyr,double& surv_ntot,nagedatlength)
// 	surv_pars = get_surv_pars(surv_seas,len_cl,crs_ind,towsperyr,surv_ntot,surv_nlengthed,qs,surv_natlength,surv_nagedatlength);
  RETURN_ARRAYS_INCREMENT();
	/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	//Returns vector "psample": Fs1, plengthed1, plength1, pagedatlength1, Fs2, plengthed2, plength2, pagedatlength2
	//
	//Requires: current survey season (surv_seas), current length class (len_cl), current cruise indicator (crs_ind),
	//current towsperyr (towsperyr), current survey ntotal (surv_ntot), current survey nlengthed (surv_nlengthed), survey qs (qs), 
	//numbers at length matrix (natlength), numbers aged at length matrix (nagedatlength)
	//
	//For now: we will condition on mles of probabilites of entering length sample, of a given length and being aged at a given length
	// as well as mles for LVB growth parameters and length-weight relationship
	/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	dvariable nlength = 0.0; 
	dvar_vector psample(1,8);
  psample.initialize();
	double naged = 0.0;
		
	if(len_cl > 0 & surv_seas > 0){
		naged = nagedatlength(crs_ind, len_cl);
		if(surv_seas == 1){				
			//spring survey
			psample(1) = qs(1) * towsperyr;
			if(naged > 0.0){
				nlength = natlength(crs_ind, len_cl);
				psample(2) = surv_nlengthed/surv_ntot;
				psample(3) = nlength/surv_nlengthed;
				psample(4) = naged/nlength;
			}
		}
		if(surv_seas == 2){
			//fall survey
			psample(5) = qs(2) * towsperyr;
			if(naged > 0.0){
				nlength = survey_natlength(crs_ind, len_cl);
				psample(6) = surv_nlengthed/surv_ntot;
				psample(7) = nlength/surv_nlengthed;
				psample(8) = naged/nlength;
			}
		}
	}
	RETURN_ARRAYS_DECREMENT();
	return psample;

FUNCTION dvar_vector get_cat_pars(int& yr, int& len_cl, dvariable& N_hat, double& cat_nlengthed, dvariable& F_cat,dmatrix& natlength, dmatrix& nagedatlength)
  RETURN_ARRAYS_INCREMENT();
	/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	//Returns vector "psample": F_catch, plengthed_catch, plength_catch, pagedatlength_catch
	//
	//Prior to use set: current year (yr), current length class (len_cl), current catch Nhat (N_hat),
	//current catch nlengthed (cat_nlengthed), current commercial F (F_cat), numbers at length matrix (natlength),
	//numbers aged at length matrix (nagedatlength)
	//
	//For now: we will condition on mles of probabilites of entering length sample, of a given length and being aged at a given length
	// as well as mles for LVB growth parameters and length-weight relationship
	/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	dvariable nlength = 0.0; 
	dvar_vector psample(1,4);
	psample.fill_seqadd(0.0,0.0);
	dvariable naged = 0.0;
	
	if(len_cl > 0){
		naged = nagedatlength(yr, len_cl);
		psample(1) = F_cat;
		if(naged > 0.0){
			nlength = natlength(yr, len_cl);
			psample(2) = cat_nlengthed/N_hat;
			psample(3) = nlength/cat_nlengthed;
			psample(4) = naged/nlength;
		}
	}
	RETURN_ARRAYS_DECREMENT();
	return psample;

FUNCTION dvariable get_s1_selectivity_age(int selectivity_type, dvar_vector& selectivity_pars, double& age)
  RETURN_ARRAYS_INCREMENT();
	dvariable selectivity_age = 0.0;
	if(selectivity_type == 0) {
//		if(age< 1.0) selectivity_age = 0.0; //assuming no fish can be caught before age 1
//		if(age>=0.0 & age < 6.0) selectivity_age = 1.0/(1.0 + mfexp(-1.0 * selectivity_pars(1)));
//		if(age>=4.0 & age < 6.0) selectivity_age = 1.0/(1.0 + mfexp(-1.0 * selectivity_pars(2)));
//		if(age>=6.0) selectivity_age = 1.0; //assuming full selectivity for 6+
		selectivity_age = 1.0; //assuming full selectivity for 6+
	}
	if(selectivity_type == 1) selectivity_age = 1.0/(1.0 + mfexp(-1.0 * mfexp(selectivity_pars(1))*(age - mfexp(selectivity_pars(2)))));
	
	RETURN_ARRAYS_DECREMENT();
	return selectivity_age;

FUNCTION dvariable get_s2_selectivity_age(int selectivity_type, dvar_vector& selectivity_pars, double& age)
  RETURN_ARRAYS_INCREMENT();
	dvariable selectivity_age = 0.0;
	if(selectivity_type == 0) {
//		if(age< 1.0) selectivity_age = 0.0; //assuming no fish can be caught before age 1
//		if(age>=0.0 & age < 2.0) selectivity_age = 1.0/(1.0 + mfexp(-1.0 * selectivity_pars(1)));
//		if(age>=2.0 & age < 6.0) selectivity_age = 1.0; // /(1.0 + mfexp(-1.0 * selectivity_pars(2)));
		selectivity_age = 1.0; //assuming full selectivity for 6+
//		if(age>=6.0) selectivity_age = 1.0; //assuming full selectivity for 6+
	}
	if(selectivity_type == 1) selectivity_age = 1.0/(1.0 + mfexp(-1.0 * mfexp(selectivity_pars(1))*(age - mfexp(selectivity_pars(2)))));
	
	RETURN_ARRAYS_DECREMENT();
	return selectivity_age;

FUNCTION dvariable get_catch_selectivity_age(int selectivity_type, dvar_vector& selectivity_pars, double& age)
  RETURN_ARRAYS_INCREMENT();
	dvariable selectivity_age = 0.0;
	if(selectivity_type == 0) {
//		if(age< 1.0) selectivity_age = 0.0; //assuming no fish can be caught before age 1
//		if(age>=0.0 & age < 6.0) selectivity_age = 1.0/(1.0 + mfexp(-1.0 * selectivity_pars(1)));
//		if(age>=6.0) selectivity_age = 1.0; //assuming full selectivity for 6+
		selectivity_age = 1.0; //assuming full selectivity for 6+
	}
	if(selectivity_type == 1) selectivity_age = 1.0/(1.0 + mfexp(-1.0 * mfexp(selectivity_pars(1))*(age - mfexp(selectivity_pars(2)))));
	RETURN_ARRAYS_DECREMENT();
	return selectivity_age;

FUNCTION void get_nhat_cohort(int& cohort_ind ) // Just set nhats (a global variable)
	/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	//Returns variable "ll_cohort"
		
	//periodinfo_ddata matrix rows: period_landmtyr, periodinfo_catch_nlengthed, periodinfo_survey_nlengthed, periodinfo_survey_ntotal,
	//periodinfo_towsperyr, periodinfo_startj, periodinfo_timeyr, periodinfo_yearlengths
	//
	//periodinfo_idata matrix rows: periodinfo_year, periodinfo_month, periodinfo_survey_season, periodinfo_cruise_indicator, 
	//periodinfo_Findicators
	
	//NOTE!!!: catch_length, survey_jcatchdate, survey_length must be in chronological (catch date) order
	//NOTE!!!: cohort_start_ind needs to be the interval where the cohort birthdate occurs!!!!
	
	nhats(cohort_ind).initialize();
	dvariable mean_wt_hat = 0.0, catch_Nhat = 0.0, F_catch = 0.0;
	dvar_vector period_probs(1,4);
	dvar_vector surv_pars(1,8);
	dvar_vector cat_pars(1,4);
	int n_catch, n_survey, current_length_class = 0;
	dvariable cum_p_surv = 1.0, current_length = 0.0, cum_page_s1 = 0.0, cum_page_s2 = 0.0, cum_page_catch = 0.0;
	double time = 0.0, current_age = 0.0;
	double current_fish_length = 0.0, current_fish_jcdate = 0.0;
	//starting index for vector of fish indices (length is # of all fish in all cohorts)
	int current_fish_ind_catch = catch_cohort_data_starts(cohort_ind);
	int current_fish_ind_survey = survey_cohort_data_starts(cohort_ind);
		
	for(int i=cohort_start_periods(cohort_ind);i<=nperiods;i++){
		//period specific info 5 idata and 8 ddata
		period_idata = column(periodinfo_idata, i);
		period_ddata = column(periodinfo_ddata, i);
		if(i==cohort_start_periods(cohort_ind)){
			if(cohort_ages_at_start_periods(cohort_ind) < 0){
				current_age = 0.0;
				//need to shorten the first interval for the birthdate of the cohort
				period_ddata(7) = period_ddata(7) + (period_ddata(6) - cohort_jbdates(cohort_ind))/period_ddata(8);
				//then change the start date
				period_ddata(6) = cohort_jbdates(cohort_ind);
			}
			else current_age = cohort_ages_at_start_periods(cohort_ind); //cohorts starting before first year of model
		}
		F_catch = F(period_idata(5));

		//number of fish caught in current period for this cohort
		n_catch = periodinfo_n_catch_aged(cohort_ind,i);
		n_survey = periodinfo_n_survey_aged(cohort_ind,i);		
		
		mean_wt_hat = sum(elem_prod(wt_at_len_hat,extract_row(catch_patlength,periodinfo_idata(1,i))));
		catch_Nhat = 1000.0 * periodinfo_ddata(1,i)/mean_wt_hat;
					
		if(n_catch > 0){			
			//get lengths to use for current period
			for(int j = 1; j <= n_catch; j++)
      {
				current_fish_length = catch_catchlength(current_fish_ind_catch);
				nhats(cohort_ind,1) += N0(cohort_ind) * get_nhat_catch_fish(period_ddata(7),F_catch,cum_p_surv,
                   current_fish_length, catch_Nhat,catch_natlength,catch_nagedatlength,
					survey_natlength,survey_nagedatlength, current_age);
				
				current_fish_ind_catch += 1;
			}
		}
		if(n_survey> 0){
			//get lengths and catchdates to use for current period
			for(int j = 1; j <= n_survey; j++){
				current_fish_length = survey_catchlength(current_fish_ind_survey);
				current_fish_jcdate = survey_jcatchdate(current_fish_ind_survey);
			
				nhats(cohort_ind,2) += N0(cohort_ind) * get_nhat_survey_fish(F_catch, cum_p_surv, 
                                 current_fish_length, current_fish_jcdate, 
					                     catch_Nhat,catch_natlength,catch_nagedatlength,
					                     survey_natlength,survey_nagedatlength, current_age);
				current_fish_ind_survey += 1.;
			}
		}
		
		current_length = Linf * (1.0 - mfexp(-k_LVB * (current_age - a0)));		
		if(current_length >= 0.5){
			for(int j = 1; j <= nlengths; j++){
				if((current_length >= double(lengths(j)) -0.5) & (current_length < double(lengths(j)) + 0.5)){
					current_length_class = lengths(j);
					break;
				}
			}
		}
		
		//not first period
		time = period_ddata(7);
		
		cat_pars = get_cat_pars(period_idata(1),current_length_class,catch_Nhat,period_ddata(2),F_catch,
			catch_natlength,catch_nagedatlength);
		surv_pars = get_surv_pars(period_idata(3),current_length_class,period_idata(4),
			period_ddata(5), period_ddata(4), period_ddata(3), qs, survey_natlength, survey_nagedatlength);
		int survey1_selectivity_type = 0, survey2_selectivity_type = 0, catch_selectivity_type = 0;
		surv_pars(1) = surv_pars(1) * get_s1_selectivity_age(survey1_selectivity_type, survey1_selectivity_pars, current_age);
		surv_pars(5) = surv_pars(5) * get_s2_selectivity_age(survey2_selectivity_type, survey2_selectivity_pars, current_age);
		cat_pars(1) = cat_pars(1) * get_catch_selectivity_age(catch_selectivity_type, catch_selectivity_pars, current_age);
		period_probs = get_probs(time,surv_pars,cat_pars);
			
		cum_page_s1 += cum_p_surv * period_probs(2);
		cum_page_s2 += cum_p_surv * period_probs(3);
		
		cum_page_catch += cum_p_surv * period_probs(4);
		cum_p_surv *= period_probs(1);
		current_age += time;
	}
	
		// cout << "nhat_catch: " << nhats(1) << endl;
		// cout << "n_catch: " << n_catch_by_cohort(cohort_ind) << endl;
		// cout << "nhat_survey: " << nhats(2) << endl;
		// cout << "n_survey: " << n_survey_by_cohort(cohort_ind) << endl;

FUNCTION dvariable get_nhat_catch_fish(double& time, dvariable& F_catch, dvariable& cum_p_surv, 
         double& catchlength, dvariable& N_hat, dmatrix& cat_natlength,dmatrix& cat_nagedatlength, 
         dmatrix& surv_natlength,dmatrix& surv_nagedatlength, double& current_age)
  RETURN_ARRAYS_INCREMENT();
	/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	//Returns variable "ll_catch_fish"
	//
	//period_ddata: landmtyr, catch_nlengthed, survey_nlengthed, survey_ntotal, towsperyr, startj, timeyr, yearlengths
	//period_idata matrix rows: year, month, survey_season, cruise_indicator, Findicators
	//
	//Requires: duration of period (time), survey q's (qs), natural mortality (m), commericial F (F_catch),
	//cumulative survival probabiltiy up to beginning of period (cum_p_surv), 
	//
	//lengths of fish at capture (cat_catchlengths), number of length classes (nlengths),
	//length classes (lengths), current catch Nhat (N_hat)
	//numbers at length matrix for catch (cat_natlength),numbers aged at length for catch (cat_nagedatlength),
	//numbers at length matrix for surveys (surv_natlength),
	//numbers aged at length matrix for surveys (surv_nagedatlength)
	/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

	dvariable nhat_catch_fish = 0.0;
	dvar_vector surv_pars(1,8);
	dvar_vector cat_pars(1,4);
	dvar_vector period_probs(1,4);
	double len, cat_nlengthed, surv_ntot, surv_nlengthed, towsperyr;
	int len_cl, yr, surv_seas, crs_ind;
	
	cat_nlengthed = period_ddata(2);
	surv_nlengthed = period_ddata(3);
	surv_ntot = period_ddata(4);
	towsperyr = period_ddata(5);
	
	yr = period_idata(1);
	surv_seas = period_idata(3);
	crs_ind = period_idata(4);;
	len = catchlength;
	for(int j = 1; j <= nlengths; j++){
		if((len >= double(lengths(j)) - 0.5) & (len < double(lengths(j)) + 0.5)){
			len_cl = lengths(j);
			break;
		}
	}
	surv_pars = get_surv_pars(surv_seas,len_cl,crs_ind,towsperyr,surv_ntot,surv_nlengthed,qs,surv_natlength,surv_nagedatlength);
	cat_pars = get_cat_pars(yr,len_cl,N_hat,cat_nlengthed,F_catch,cat_natlength,cat_nagedatlength);
	int survey1_selectivity_type = 0, survey2_selectivity_type = 0, catch_selectivity_type = 0;
	surv_pars(1) = surv_pars(1) * get_s1_selectivity_age(survey1_selectivity_type, survey1_selectivity_pars, current_age);
	surv_pars(5) = surv_pars(5) * get_s2_selectivity_age(survey2_selectivity_type, survey2_selectivity_pars, current_age);
	cat_pars(1) = cat_pars(1) * get_catch_selectivity_age(catch_selectivity_type, catch_selectivity_pars, current_age);
	period_probs = get_probs(time,surv_pars,cat_pars);
	
	nhat_catch_fish = cum_p_surv * period_probs(4);
	RETURN_ARRAYS_DECREMENT();
	return nhat_catch_fish;

	
FUNCTION dvariable get_nhat_survey_fish(dvariable& F_catch, dvariable& cum_p_surv,double& catchlength,
         double& jcatchdate, dvariable& N_hat, dmatrix& cat_natlength,dmatrix& cat_nagedatlength,
         dmatrix& surv_natlength,dmatrix& surv_nagedatlength, double& current_age)
  RETURN_ARRAYS_INCREMENT();
	/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	//Returns variable "ll_survey_fish"
	//
	//period_ddata: landmtyr, catch_nlengthed, survey_nlengthed, survey_ntotal, towsperyr, startj, timeyr, yearlengths
	//period_idata matrix rows: year, month, survey_season, cruise_indicator, Findicators
	//
	//Requires: natural mortality (M), 
	//cumulative survival probabiltiy up to beginning of period (cum_p_surv), 
	//length at capture (catchlength), catch date (jcatchdate),
	//number of length classes (nlengths), length classes (lengths), start of period (startj), length of current year (yrlen), current year (yr), 
	//current catch Nhat (N_hat), number lengthed in catch (cat_nlengthed) , commericial F (F_catch),
	//numbers at length matrix for catch (cat_natlength),numbers aged at length for catch (cat_nagedatlength),
	//current survey (surv_seas), current cruise (crs_ind), towsperyr, total number caught in survey (surv_ntot), 
	//number lengthed in survey (surv_nlengthed), survey q's (qs), numers at length matrix for surveys (surv_natlength),
	//numbers aged at length matrix for surveys (surv_nagedatlength)
	/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

	dvariable nhat_survey_fish = 0.0;
	dvar_vector surv_pars(1,8);
	dvar_vector cat_pars(1,4);
	dvar_vector period_probs(1,4);
	int len_cl, yr, surv_seas, crs_ind;
	double len, cat_nlengthed, surv_ntot, surv_nlengthed, towsperyr, startj, yrlen, time;
	
	cat_nlengthed = period_ddata(2);
	surv_nlengthed = period_ddata(3);
	surv_ntot = period_ddata(4);
	towsperyr = period_ddata(5);
	startj = period_ddata(6);
	yrlen = period_ddata(8);
	
	yr = period_idata(1);
	surv_seas = period_idata(3);
	crs_ind = period_idata(4);;
	
	
	
	len = catchlength;
	for(int j = 1; j <= nlengths; j++){
		if((len >= double(lengths(j)) - 0.5) & (len < double(lengths(j)) + 0.5)){
			len_cl = lengths(j);
			break;
		}
	}
			
	time = (jcatchdate - startj)/yrlen;
			
	surv_pars = get_surv_pars(surv_seas,len_cl,crs_ind,towsperyr,surv_ntot,surv_nlengthed,qs,surv_natlength,surv_nagedatlength);
	cat_pars = get_cat_pars(yr,len_cl,N_hat,cat_nlengthed,F_catch,cat_natlength,cat_nagedatlength);
	int survey1_selectivity_type = 0, survey2_selectivity_type = 0, catch_selectivity_type = 0;
	surv_pars(1) = surv_pars(1) * get_s1_selectivity_age(survey1_selectivity_type, survey1_selectivity_pars, current_age);
	surv_pars(5) = surv_pars(5) * get_s2_selectivity_age(survey2_selectivity_type, survey2_selectivity_pars, current_age);
	cat_pars(1) = cat_pars(1) * get_catch_selectivity_age(catch_selectivity_type, catch_selectivity_pars, current_age);
	period_probs = get_probs(time,surv_pars,cat_pars);
			
	nhat_survey_fish = cum_p_surv * period_probs(1);
	if(surv_seas == 1){
		nhat_survey_fish *= surv_pars(1) * surv_pars(2) * surv_pars(3) * surv_pars(4);
	}
	if(surv_seas == 2){
		nhat_survey_fish *= surv_pars(5) * surv_pars(6) * surv_pars(7) * surv_pars(8);
	}
	RETURN_ARRAYS_DECREMENT();
	return nhat_survey_fish;
	
REPORT_SECTION
  //OUTPUT RESULTS TO FILE "fsctred.REP"

  report << "Redfish finite-state continuous-time age-structured model" << endl;
  report << "OBJECTIVE FUNCTION VALUE: " << total_ll << endl;
  report << "N0: " << endl;
  report << N0 << endl;
  report << "qs: : " << qs << endl;
  report << "M: " << M << endl;
  report << "F: " << endl;
  report << F << endl;
  report << "survey1_selectivity parameters: " << endl;
  report << 1.0/(1.0 + mfexp(-1.0*survey1_selectivity_pars)) << endl;
  report << "survey2_selectivity parameters: " << endl;
  report << 1.0/(1.0 + mfexp(-1.0*survey2_selectivity_pars)) << endl;
  report << "catch_selectivity parameters: " << endl;
  report << 1.0/(1.0 + mfexp(-1.0*catch_selectivity_pars)) << endl;
  report << "first ll for survey: " << ll_first << endl;
  report << "misc data for first survey: " << temp_ddat << endl;
  report << "expected numbers caught in fishery and survey: " << endl;
  report << nhats << endl;
  //END OF MODEL
